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We investigate two-fluid BGK kinetic methods for binary fluids. The developed theory works for 

asymmetric as well as symmetric systems. For symmetric systems it recovers Sirovich's theory and 

is summarized in models A and B. For asymmetric systems it contributes models C, D and E which 

are especially useful when the total masses and/or local temperatures of the two components are 

greatly different. The kinetic models are discretized based on an octagonal discrete velocity model. 

The discrete-velocity kinetic models and the continuous ones are required to describe the same 

^^' hydrodynamic equations. The combination of a discrete-velocity kinetic model and an appropriate 

^ Cy , finite-difference scheme composes a finite-difference lattice Boltzmann method. The validity of the 

^>. ■ formulated methods is verified by investigating (i) uniform relaxation processes, (ii) isothermal 

' ' Couette flow, and (iii) diffusion behavior. 
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I. INTRODUCTION 



Gas kinetic theory plays a fundamental role in understanding many complex processes. To make solutions possible, 

many of the kinetic models for gases are based on the linearized Boltzmann equation, especially based on the BGK 

approxiniation'15. Even thus, only in very limited cases are analytic solutions available. Basically speaking, there are 

two options to simulate Boltzmann equation systems. First, one can design procedures based on the fundamental 

^ , properties of rarefied gas alone, like free flow, the mean free path, and collision frequency. Such a scheme does not 

Q ■ need an a priori relationship with the Boltzmann equation, but the scheme itself will reflect many ideas and/or 

O [ concepts used in the derivation of Boltzmann equation. In the best case, such a simulation will produce results being 

consistent with the solution of Boltzmann equation. The second option is to start from the Boltzmann equation and 

design numerical schemes as accurate as possible^. The discrete Boltzmann equation approach or lattice Boltzmann 

^ : method (LBM) has been becoming a viable and promising scheme for simulating fluid flowsjSJ. 

CNj LBMs for single-component fluids have been well studied, while for binary mixtures still need more clarification Ul. 

For binary fluids , although various LBMs have been proposed_p, Il[lllla[l3,[lll[ll[ll[ll[ll ElEl 111 111 llfl 

miH , most of them |1 11 111 III 111 111 111 111 111 [llM 113 are based on the single-fluid theoryl2g|~^For systems 

with different component properties, a two-fluid theory is necessary. Sirovich's two- fluid kinetic theory '24] works for 

^i ■ (approximately) symmetric systems where the two components have (approximately) the same total masses and local 

(^ I temperatures. A LBM based on Sirovich's theory and for the complete two-dimensional Navier-Stokes equations(NSE) 

"p-; , is given in ^Z|. This LBM is based on a two-dimenaional model with sixty-one discrete velocities (D2V61). Many 

C^ compressible fluids can be well described by the Euler equations ,25]. In fluid mechamics of low-speed flow, the 

G . temperature remains nearly constant and consequently the isothermal NSE description is extensively used|25|. From 

I ' the Chapman-Enskog procedure |2a| the Euler equation is a lower-order approximation compared with the NSE. The 

' O isothermal NSE is a simplifled case of the complete NSE. For the above two kinds of systems, using the LBM for 

complete NSE system is not neccessary and computationally inefficient. In this study we generalize Sirovich's theory 

so that it works also for asymmetric systems where the total masses and/or local temperatures of the two components 

are greatly different, then formulate LBMs for the two kinds of systems. The LBMs formulated here require simpler 

discrete velocity models(DVMs). For the Euler-equation system a DVM with thirty-three discrete velocities (D2V33) 

is enough. For the isothermal NSE system, a D2V25 is sufiicient. 

^ This paper is arranged in the following way: In section II we review and develop the two-fluid BGK kinetic theory. 

5^ , Sirovich's original treatments are clarifled and summarized in models A and B. For asymmetric systems three kinetic 

models (C, D and E) are derived. The hydrodynamics and diffusion behavior of the model systems are discussed. 

In section III the kinetic models are discretized based on a multispeed discrete velocity model. Then, possible FD 

schemes are given and the corresponding numerical viscosities and diffusivities are analyzed. Numerical tests are 

shown in Section IV. Section V concludes the present paper. 

II. TWO-FLUID BGK KINETIC THEORY 

In a binary system with two components, A and B, roughly speaking, the approach to equilibrium can be divided 
into two processes. One is referred to as Maxwellization (i.e., each species equilibrates within itself so that the local 



distribution function approaches to its local Maxwellian) . The other is the equilibration of species (i.e., the differences 
in hydrodynamic velocities and local temperatures of the two components eventually vanish) . Correspondingly, the 
interparticle collisions fall into two categories: self-collisions (collisions within the same species) and cross-collisions 
(collisions between different species') |5ll24). 



A. General description 



For a two-dimensional binary gas system the BGK kinetic equations readQ, 



dtr + ^""—f + a^ • Tfr/-^ = J-^-^ + J-^^, 



dr 
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f^ ( /^) s-iid v"^ (v'^) are the distribution function and particle velocity of the component A {B); /-^C^) and J'4B(o) g^j.g 
the local Maxwellians which work as references for the self- and cross-collisions; n^, u^, T^ are the local number 
density, hydrodynamic velocity and temperature of the species A; u, T are the local hydrodynamic velocity and 
temperature of the mixture after equilibration process; a^ is the acceleration of the species A due to the effective 
external field. For species A, we have 



n- = dv-r , 



(7) 
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where p and e|7it (^o ) £^re the local mass density and internal mean kinetic energy (hydrostatic pressure) of species 
A, ks is the Boltzmann constant. For species B, we have similar relations. 
For the mixture, we have 



n = n^ + n^, p = p-^ + p^, 



(12) 
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(13) 



nksT 



dv^l (v-^ - ufm^f^ + fdv^\{.^ ufm^f^, 



(14) 



eint(^o) = nksT, 



(15) 



where n, p, u, T, eint, ^o Sive the total number density, total mass density, barycentric velocity, mean temperature, 
total internal energy, and total hydrostatic pressure, respectively. It is easy to find the following relations, 
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Here three sets of hydrodynamic quantities [{p^, u^, T^), (p-^, u"^, T^) and (p, u, T)] are involved. If assume that 
the two components are in local equilibrium, implying that T^, T^ , T can be replaced by T^°) and u"^, u-^, u can be 
replaced by u^"^ in the definitions of /'^(°) , /-^(o) and f^B(o)^ ^^ arrive at the one-fluid theory and Eq. (|17|l recovers 
Dalton's law'25'|, where T^^^ and u^°^ are the temperature and the velocity of the system in the complete equilibrium. 
It is clear that the one- fluid theory is conditionally valid. If the differences among T'^, T^, T and/or among u'^, u^, 
u are not small, the above replacements result in large errors. Since each set of the hydrodynamic quantities can be 
described by the other two sets, in such cases, a two-fluid theory is preferable. Without loss of generality, we require 
the description to be dependent on {p^, u^, T^) and (p^, u^, T^)'2i\. 

A key point to complete the two-fluid kinetic description is how to calculate the local Maxwellian /^^(o) (/^"^(°)). 
Within Sirovich's original treatments, it is Taylor expanded around /"*(") (/^''°'') to the first order of flow velocity and 
temperature[2J|. This treatment is reasonable when the hydrodynamic properties of the two components are nearly 
symmetric, i.e., lu — u'^l « lu — u^l, \T — T"*! « IT — T^l. To make a general theory working also for asymmetric 
systems where the hydrodynamic properties of the two components are greatly different, we introduce the reference 
distribution function in a general way and do the Taylor expansion around it. 
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where the second superscript "r" means "reference" and the corresponding quantities are the reference hydrodynamic 
quantities which take values in the following way. 
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Let us make the solutions more explicit. Firstly for u, from Ea. l|13(l we have 



u — u 



n^m^ 



Then for T, from Eq. ^, when T"^ > T^ we have 
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when T^ < T^ we have 
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Considering together ipUjI - iPH) gives 
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In the case of u^'' = u^ and T^"" = T^, 5-^^°) gets back to /^(°l 

Both of p^ and T^ are local quantities. Their values are functions of position and time. It is possible for such a 
phenomenon, p^{ri,t) > p^{ri,t) but p^{r2,t) < p^{r2,t), to occur, where ri and r2 are two different positions in 
the system. While in a theory it is not convenient to use the reference state in such a way: u"^'"(ri, i) = u'^(ri, t) and 
u'^''(r2, i) = u^(r2, t). Instead, we prefer to use one of the two possibilities, u'^''(r, t) = u"^(r, t) or u"^'"(r, t) — u^(r, t), 
in the whole system, where r is an arbitrary position in the system. For T"^'' we have the same preference. To that 
aim, p"^,p^ ,T^ and T^ in the criteria (|25|) and (|26|1 are replaced by their spacially averaged values, p^, p^ , f^ and 
T^ , respectively. This treatment is reasonable from a statistical sense. 

B. Kinetic models for symmetric systems 

For systems with p^ « p^ and f^ fa T^ , we can use g^'-°'' = f^'-°\ g^'-°'' = f^^°\ i.e, Sirovich's kinetic theory. In 
this case the equations for the two components are symmetric. The cross-collision term in (^ becomes 
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where 
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If wc concern the hydrodynamics only up to the NSE level, f^ (/^) in the force term can be replaced by /"^("^ (Z^*-"^). 
The BGK model (I1I6|I can be rewritten as 
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where 
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the expressions of Q^^ and Q^"^ are obtained from Eqs. (|31|l and H32|) via formal replacements of the superscripts A 
and B. In the isothermal case, T"* = T^ = T, the expression of Q^^ is simplified as 



Q^'' = -^^^i{-''-n^)■in^-n^)■ 



(33) 



For the convenience of description, the kinetic model with H29|) -()32 |) is referred to as kinetic model A; the one with 
(1231- (ED: (ESI is referred to as kinetic model B. 



C. Kinetic models for asymmetric systems 

1. Kinetic model C: for isothermal systems with p^ > p^ 
For such a system, g-^C^) — f'^^°\ and 
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2. Kinetic model D: for systems with p > p^ and T > T^ 



The references are g^^^^ = /^(o) and 
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where 
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3. Kinetic model E: for systems with p > p and T < T 
In this case, the reference velocity and reference temperature for both /^^C) and fBA{o) ^^^^ ^A ^^-^^ yS^ respectively. 
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Within the kinetic model E 



S(0) ^ P 



jB ^^B _ ^A 
2kBTB 



21 



Q 



)AA 



BB 



rAA 



rBB 



jA _ fA(0) 
jB _ yS(0) 



1 


^AB 


1 r 


tBA [ 



J A _ gA(a) 

fB^gBiO) 



(49) 

(50) 
(51) 



Q 



AB 9 



A(0) 



e^r 



,-(v--u-).(u--u-) + -^ 



T^^nm^ 



A ,,A\ 



(v^ - u 
29^ 



B rpA\ 



{T" -T 



~M' 



29^^ 



A ,,B\2 



{W^ - U 



(52) 



Q^^ = 






(v^ 



29^ 



- 1 



(T^ - T-^) 



-M'' 



(v^ 



u 



26 



«s 



(u^ 



,A\2 



(53) 



D. Hydrodynamics and diffusion 

1 . Hydrodynamics 



A connection between a kinetic model and corresponding hydrodynamics is the Chapman-Enskog analvsis|26J. All 
above kinetic models contribute to (i) the same continuity equation at the Euler and the NSE levels, 






(ii) the same Euler momentum equations 
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and (iii) the same NSE momentum equation for component A^ 

d dPf^ 



d_ 
dt 



^ (/'""") +£:(A^^') + ^^"' "'""''" 



AB 



dr/j 



drp 



P < + 



AB 



where 



AB 






(57) 



(58) 



(59) 



describes the momentum transferred from component A to B, and it is also the diffusion flux density which will be 
clear from a later equation H72|) : 



is the stress tensor, 
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is the viscous stress tensor, and 
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is the viscosity. 

Moldes A and B contributes to symmetric hydrodynamics for the two components. The Euler energy equation of 
model A for component A reads. 
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where 
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is the local total energy, and 
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is the heat transfcred from component A to _B. 

The NSE momentum equation for component B from model C reads 
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where the definition of P^a is similar to that of P^a, and 
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is an additional stress tensor due to the asymmetry of densities of the two components. 

The Euler energy equation of model D for component A is the same as Eq. (|63ll and for component B reads 
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where the definition of e^ is similar to that of e and 
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The Euler energy equations from kinetic model E are as follows, 
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2. Diffusion 



From the continuity equations (|54|I - H55|I we have 
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where j^^ is given in Eq. H59() and it is the amount of the component A transported relative to the component B by 
diffusion through unit area in unit time. For the incompressible fluids where p is a constant, the continuity equation 
(f7^ is equivalent to the following diffusion-convection equation, 

dip d , , d 



where Lp = p / p. The diffusion velocity 
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^) is determined by the monientum equation. We can find a simple 



dt ■ ar„(^"") = ~a^L^(^ 



relation for it in the following case: We consider a binary system without external forces and where the flow velocities 
u^, u^ are small and their derivatives can be regarded as higher-order small quantities. From the momentum Eq. 
(|56|l or (|58|l . by neglecting the second and higher-order terms in u^ and/or wf , then using the definition (|11|) . we 
obtain 
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If further assume the system to be isothermal, the density flux of component A reads 
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is the diffusivity of component A. Eq. (|75|l is Fick's first lawj28J. 
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From Eqs. (|73|l and (|75|) we have 
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is a constant, the diffusion-convection equation l|78() 



(79) 



In the case where the barycentric velocity field is zero and t^^ 
reduces to Fick's second law^2S]. 

dp_ ^j^A 9 dip 
dt dra dra 

Under the present treatment cross-collisions contribute to the viscous behavior and are responsible for the inter- 
diffusion as well as momentum and heat exchanges between the two components. The momentum and heat exchanges 
between the two components occur not o nly at the Navier-Stokes level but also at the Euler level 29] , which is different 
from the case in the one- fluid theory[l^|23, but consistent with the two-fluid relations H16() and H17|l . 



III. DISCRETE KINETIC MODELS 



A. General description 



Based on the following discrete velocity model, 

vo = 0, Vfci = Vk 
the kinetic equations read 
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where subscript k indicates the fc-th group of particle velocities and i indicates the direction of the particle speed. 
The DVM H8U|I is isotropic up to its seventh rank tensors [3(|. The discrete kinetic model, (I81|l - (|82|l . is required to 
recover the same hydrodynamic equations as those of its continuous version. This requirement is used to formulate 
the multispeed-discrete-velocity kinetic models. 



B. Models for isothermal and compressible Navier-Stokes equations 

1. Discrete-velocity kinetic model B 
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Due to the symmetry of the two components, we show results only for the component A, 
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The Chapman-Enskog analysis shows that, to get the isothermal NSE equations, 1)54(1 and (|58|l . the following require- 
ments on the discrete equilibrium distribution function. 
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are necessary and also sufficient. 

The requirement H89|l contains the third order of the flow velocity u"^. So it is reasonable to expand fj^^ ' in the 
polynomial form to the third order in the flow velocity. 
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The left-hand side of Eq. ((89|1 with the truncated /j,/ has sixth rank tensor in particle velocity v^. Therefore, to 
recover the correct hydrodynamical equations, the based DVM should be isotropic up to its sixth rank tensor. DVM 
(|80|l satisfies the need. 

To satisfy ()86|l. we require 



E^"-i' 



(92) 



ki 



E 

ki 



FM^,vLufui = 9^ {u^) 



(93) 
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To satisfy (|87|l . we require 



Y. Fkvto.vi^uf = e 



^u^ 



ki 



Y. F^vi^vi^vi^vi^ufuiu^ = 3 (e^) (u^) 



^^u^. 



ki 



To satisfy (jHHl, we require 



YFkvicvip^Q^Sc^ 



ki 



Y.PkVk^A4^ivL4''^^i^^) K) '^"/3 + 2"f 



A„,A 



ki 



To satisfy (|89|l . we require 



Y.^k^k^o.^t^pvt^■y^t^^H = (0^)' K^Pl + "^-^^a + ^-^^/j) 



(94) 
(95) 

(96) 

(97) 

(98) 



EF^7I^ 71^ 71^ 71^ 71^ lA 7/^7;'^7;^ 
-^fc "kia"kip"ki-f"ki^"kiiT"kir^'^(, "tt "jj 



^\'^„,^„,^„,A 



3 (e^) (w-^) {ut5fi^ + w^<5-ya + u;^<5a0) + 6 (e^) u^uj 



(99) 



If further consider the isotropic properties of the discrete velocity model, the above 8 requirements reduce to the 
foUowing four ones. Requirement H92|l gives 



Requirements (jHSll, IMIl, 1® give 



Requirements (jHEl), JHHl, (EHl give 



Requirement (|^ give 



E^^" = i- 



fci 



k 



E^^"K)' = (0^ 



E^'^l-^ 



i' = 6fe^)'. 



(100) 



(101) 



(102) 



(103) 



To satisfy the above four requirements, four different particle velocities are sufficient. We choose a zero speed, Vq = 0, 
and other three nonzero ones, u^ (fc — 1, 2, 3). From (|101(I - H103() it is easy to find the following solution. 






'^t = Q 



e^{K+i 



vt2) -^Q^ 



<+i + < 



'-•k+i 



k+2) 



24 



(e-)^} 



$^ = 4 



K)' {«1<2)' - K)' [K+l)' + (-^+2)'] + (V^)"} 



where k = 1,2,3 and vf — v^, v^ — v^. From pOO|) we get 



Fo^^l-8^F, 



(104) 

(105) 
(106) 

(107) 



fe=i 
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2. Discrete-velocity kinetic model C 



The description for component A is the same as that in discrete-velocity kinetic model B. For component B, 

1 

D D ,/ hi J h'i 



nBB _ 1 \ rB _ fB(0) 

^ki ^ ^BB ■'ki Jki 



^BA 



fB _ B{0) 
■I ki fJki 



(108) 



Qki 



BA 

ki 



S(0) 
3ki B* ( B 

QB ^-D V'ki 



(u^ - U^). 



(109) 



/j.j and /j.j are formulated in the same as those in discrete model B. Additionally, g^.^ should be formulated in 
a similar way. Due to similar reasons, gj^- is expanded as 



s(o) 
9ki 



Br'B 



n"Gt 



A\2 



29^ 



QB 



A\2 



1- 






vfipuf + 



_,,-B .,B A A 



6(6^) 



'S'^kii^kiir'^kiri'^i "tt "">) f + 



(110) 



where 



^ 2^QB 



exp 



(-f) 



29^ 



(111) 



The formulae for G^ can be obtained through formal replacements in Eas. p04|l - H107|) : 9"^ —* Q^ = ksT /m^ ^ F^ ~> 

C. Models for the complete Euler equations 

LBMs for s ing le-component Euler equation have been constructed by several authors. [See Yan, et al 31] and 
Kataoka, et al|32| for examples.] In this section we formulate the discrete- velocity kinetic models A, D and E for the 
complete Euler equations of binary fluids. 



1. Discrete-velocity kinetic model A 
The equations for component A are the same as Eqs. H83|l . (|85|l with 



QtP - 



^aT 1 ^^ V 



A (,,A 
ki 



U 



-M' 



29-4 



(u^-u^)+M^ 



A ,.B\2 



K»-u^^' 



29^ 



- 1 



(T^ „ T^) 



(u^ - u 



(112) 



The Chapman-Enskog analysis ^6j shows that, to recover the same Euler equations, 154|l . (|56l) and (|63|l . besides 
(|SS)l - (|S^ . one more requirement 



Eo'^'^K)'^^-^^*/^/- 



A „,A fA{0) 
ki 



ki 



2P^Q^5^p + ^ {u^fs^^ + iP^utuj + \p^ {u^'fut 



(113) 



is necessary. Correspondingly, f^^ should be Taylor expanded to the fourth order of flow velocity, the DVM should 
be isotropic up to its seventh rank tensor. Again, DVM (|8U|1 satisfies the need. To satisfy (|113f) . we require 

\2 



EuA \^k ) „,A ^,A _9/'oA^^r 
^k — ^kia^kip - 2 (9 ) dap , 

ki 



(114) 
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J2Pk^^^^^o.v^^P<^^^^''t< = 3 {e^f [{u^)'s^p + 2uiui 



ki 



(115) 



ETT'A y^k ) A A A A A A A A A A 
^k ^ ^kia'^kifi'^kii'^kiTT'^kiri^kiX'^e, '^n'^rj'^X 



= 12 (e-4)' [{u^ySo^p + 4 {u^)\y^ . (116) 

Finally, we have five requirements on Fj^. Four are shown in Eqs. pUU|l - pU3|) and the fifth is 

J2f,^{v^)'^48{Q^)\ (117) 

k 

To satisfy the above five requirements, five particle velocities are sufficient. We choose a zero speed, Vq — 0, and 
other four nonzero ones, v^ {k = 1, 2, 3, 4). It is easy to find the following solution, 



F„^ = 1-8^F,^, 



fc=i 



F' 






where 



3 3 






HvtfnuUv^f-ivU" 



k = 1,2,3,4, and v^_^_, = v^, (j = 1, 2, 3). For component B we have similar results. 



(118) 
(119) 

(120) 
(121) 



2. Discrete-velocity kinetic model D 



The equations for component A arc the same as those of model A and for component B are as follows, 



Q'' = -^ 



Jki J ki 



S(0) 



1 



^BA 



fB B{Q) 

Jki 9kt 



(122) 



Qf.^ - 



5(0) 

9ki 

QBr 



-PF 



/ig* (vf, - u^) . (u^ - u^ 



Mt 



26^'- 



-1 



(rpA _ rpl 



(v 



ki 



29^'- 



(u^ - u^f 



(123) 



/j,^ and /j,^ are formulated in the same way as in model A. g^,^ is expanded in the similar way to (|110|) but to 
the fourth order of flow velocity. The formulae for G^ can be obtained through formal replacements in Eas. (|118|) - (|121|l : 
e^ ^ QBr ^ kBT^/m^, F^ -^ Gf . 



3. Discrete-velocity kinetic model E 



Within model E 



iki 



^BB 
iki 



rAA 



rBB 



fA _ f^ 

J ki J k 
Jki ~ Ik 



A{0) 



rAB 



rBA 



f^ — n 
Jki ilk 

fB 

Jki - 9k 



MO) 



(124) 
(125) 
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iki 



A(0) C 



« 



-M 



2QAr 



1 



(u 



A ^B^2 



2eAr 



-1 



(T^ - T^) 



(126) 



^BA 

Iki 



-B(O) f 

5fe2 J .,B* ( B 

QB ] ^^D [^ki 



2QB 



A\2 



(u« - u-^) 



B* 



u^ •K-u^)+Mt 



{< - u^) 



26 



«s 



^ThA _ Tni 



(127) 



rA{0) rB{0} 



A(0) 



B(0) 



fki ' /fci ^'"^ formulated in the same way as in model A. The formulations of 5^^ and g^.^ are similar to 
those in model D. The requirements on 5^^ can be obtained from those of model D by using formal replacements: 



A c,Br 



replacements; v^ 



QB . Then the requirements on 5^,^ can be obtained from those on 5^^ by using formal 



e^ 



cvAr 



e 



kBT^/m^, gl 



S(0) 



ilki ■ 



D. Finite-difference schemes, spurious viscosities and diffusivities 



The time evolution of a discrete-velocity kinetic model can be solved numerically by using appropriate finite- 
difference schemes. There are various options for calculating the time derivative and the advection termp^Tl Ha. l3^ . 
In a practical simulation the real evolution equation of /^ is not Eq. (|81|l but 



dt 



r)2 
Q^2Jk^ 



da 



kia. 



■ ^«fe»« 






}t^ + ' 



-.AB 

iki 



{< 



U^) rA{0) 



QA 



-r 

Jki 



(128) 



where smaller terms in the second and higher orders of At or Aa have been neglected; the factors 9 and '0 can be 
specified for various FD schemes. The extra terms in 9 and tp contribute to the spurious viscosities and diffusivities 
in the simulation results. To check the spurious viscosities and diffusivities, one needs do again the Chapman-Enskog 
analysis to Eq. (|128|l and compare the hydrodynamic equations with those of the continuous models. The recovered 
mass and momentum equations from H128() are 



dp^ d , A A\ 



dj^ 



AB,S 



dra 



(129) 






d 

dra 



(p'^Uo 



dra 



br+j 



AB,S1 



(130) 



and 



di ^^^"") ^d^ (p^^^p) 



d 

drf3 



-4 



drp 



dt 



■V 



dr/sj 



d 
drp 



dPl 
dra 



P^at + I 



AB 



n 



ap 



{p^utuj) 



dP,\ 



AB 



dt 



aAB 

Ja 

^AB 



P^ai 



(131) 



where 



3a 



AB,S 



= -\(^7r.-^i^]{p^<) 



dt 



dra 



(132) 



is the spurious diffusion flux density, d/drp (d/drp) = d'^/dr^. The spurious diffusivity and viscosity are coupled in 
the real momentum equation p31|) . The real momentum equations for component B and the real energy equations 
can be considered in a similar way. 
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Which FD scheme to use depends on the question under consideration. Since the higher-order schemes for time 
derivative require more memory, the forward Euler scheme is generally used. In binary systems concentration gradients 
drive the diffusion behavior. For systems with large density gradients, the space centered scheme is less stable and 
the wiggle phenomena of the second-order upwind, the Lax-Wendroff and the Beam- Warming schemes introduces 
unphysical oscillations of fluid densities[i3,[lMl23|- Therefore, for such a system, the first-order upwind scheme 



df, 



da 



Jk 



i.I Jki 



-Ac 






>o 

<0 



(133) 



is generally preferred, where the third subscripts / — 1, /, / + 1 in Eq. (|133|) indicate consecutive mesh nodes in the 
a direction and Aa is the space step. In such a FDLBM scheme, 6 = At/2, ip = —Aa/2 if w^^^ > and -0 = Aa/2 if 



<a < fl- 



it should be noted that besides the FD schemes and truncation errors of the machine, the numerical errors from 
the DVMs also contribute to spurious diffusivities and/or viscosities. The smaller the hydrodynamic velocity, the less 
this part of contribution. Other discussions on the origin of spurious velocities and possible remedies are referred to 

[milillSllSlElSElEi 



IV. NUMERICAL TESTS 

As mentioned above, in a practical simulation the numerical errors have three resources, the formulated DVM, 
the spacial discretization and the time discretization. We first check the case where the spacial FD scheme has no 
contribution - the uniform relaxation process where the physical quantities are only functions of time. For the velocity 
equilibration, the five kinetic models give the same expression. 



d{i 



at 



P 



rAB 



^BA 



(134) 



For temperature equilibration, model A gives 

d{T^-TB) 
dt 



^AB 



+ 



2kBnp 



rBA 



(rpA _ rp± 



^BA 



-AB 



K-u^)% 



(135) 



model D gives 



g (J.A _ J.B) 



dt 



1 / n^ 
n \T 

1 
2kBnp 



A B 
P P 



1 1 



rAB 



rBA 



rBA 



[P' 



P") 



(U^ 



U^)' 



(136) 



and model E gives 



d{T' 



T' 



at 



n 



1 

n \T 

1 
2kBnp 



AB 



rBA 



p^p"" 



{T' 



rAB 



rpB 

rBA 



rBA 



u-y 



(137) 



Numerical examples are shown in Fig.l. In Fig. 1 (a) we show two cases where p — p^', in the isothermal case kinetic 
models A and B are applied, while in the case of T^ « T^ only kinetic model A is applied. Fig. 1 (b) shows cases 
where p^ ^ p^ so that models A and B do not work and one has to resort on models C, D and E. For the velocity 
equilibration procedure, under the accuracy of the calculations, models A and B give the same results, models C, D 
and E give the same results. All the numerical results in (a) and (b) agree well with the theoretical ones. 

Secondly, we check a case where the advection terms make effects and viscosities exist. We use the two-fluid 
FDLBMs A and B to investigate the isothermal Couette flow for single-component fluid. The two walls, locating at 
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y = zLD/2, start to move horizontally with velocities ±t/ at i = 0, where D is the distance between the two walls. 
The simulation results of the velocity profiles agree well with the following theoretical one, 

. = 7y-f:(-lF-^^exp(-^t)sin&, (138) 

where u is the horizontal velocity, 7 = 2U / D the imposed the shear rate, j an integer . ( An example is referred to 
Fig. 2.) ^ ^ 

Thirdly, we investigate the diffusion behavior in a one-dimensional system. To make valid the relation (|77|l and 
make less the numerical errors from the spacial FD scheme, we assume that (i) the two components have equal particle 
masses 1, (ii) the initial hydrodynamic velocities of the two components are zero, (iii) the system is isothermal with 
temperature T = 1 , (iv) the initial density profiles of the two components are 

4_fl-Apifx<0 B_/1 + Apifx<0 
P ~\l + Apifa;>0'^ ^\1-Apifx>0' ^ '^•'> 

where < Ap < 1, and (v) the viscosities of the two components are small enough. Thus, the barycentric velocity 
field of this system is globally zero, V^ — r'^^, V^ — t^^, and the evolution of the density profiles follows 

"'-'^ ^'"' (vii) ■'°''- ^'"' {7m) - '"»' 

To make the numerical tests practical, when choose parameters for simulations, the following points should be con- 
sidered: (i) The accuracy of the forward Euler scheme is in the order of At and that of the upwind scheme H133|) is 
in the order of Ax; (ii) If the physical values of V^ and V^ are too small, they may be submerged by the numerical 
diffusivities. Numerical tests show that LBMs A and B can recover density profiles which agree well with Eq. (|14()|l . 
An example is shown in Fig. 3. A set of density profiles for the component A are shown in (a). To help evaluate the 
numerical errors from the DVM, the corresponding profiles of diffusion velocity u^ are shown in (b). The diffusion 
velocity u has its maximum value at x = 0. Its magnitude decreases with time. For the earliest time (t — 0.01) 
shown in this figure, |u"^| « 0.028. The numerical errors for f^,-^ ' are in the order of (m"^) for LBM B and in the 
order of {u'^f for LBM A. 

V. CONCLUSIONS AND REMARKS 

Sirovich's original two-fluid BGK kinetic theory works for symmetric systems where the two components have 
approximately the same total masses and local temperatures. This theory is clarified and generalized to describe both 
symmetric and asymmetric systems. Corresponding to different situations five kinetic models are formulated. Based 
on an octagonal discrete velocity model the five models are discretized. The discrete-velocity kinetic models and the 
continuous ones are required to recover the same Euler and/or Navier-Stokes equations. A discrete-velocity kinetic 
model and an appropriate finite-difference scheme compose a FDLBM. The formulated kinetic models work also for 
binary mixtures with disparate particle-mass components. Which model to use depends on the mean temperatures 
and the mean mass densities of the two components. 

In the present two-fluid treatment, the relaxation times of the cross-collisions contribute to both the viscous and 
diffusive effects. The interfacial tension is another aspect of thermodynamic interaction between component fluids. 
Investigating the interfacial tension is crucial in the industrial context for controlling the size and phase stability 
of mechanically dispersed droplets and other transient structures formed in the course of phase separation. For 
immiscible fluids, one way to introduce the interfacial tension is through modifying the pressure tensors [l^ by taking 
into account the the interparticlc interactions. One possibility of incorporating the interparticle interaction is through 
modifying the force terms in Boltzmann equations [l^j. In such a case, the force terms in the BGK kinetic models are 
responsible for the phase separation and interfacial tension. The acceleration a'^ is determined by the interparticle 
interactions and the external fleld. The determination of the specific form of a depends on the system under 
consideration. An interesting point is that the incorporation of the force term in the Boltzmann equation makes no 
additional requirement on the formulation process of the FDLBM. So the specific forces can be directly considered 
under the same frame. A different attempt to introduce the interfacial tension is to start from the Enskog equations 
for dense gases |^. 
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FIG. 1: Velocity and temperature equalibrations in uniform relaxation processes. The common parameters used in (a) are 

n^ = l,n^ = 2, m-^ = 2, m^ = 1, r""^ = r^^^ 

and T = 1 for the isothermal case; T"^ = 1.2, T^ = 0.8, u-*'"' = u^'"' 



— 1. Additionally, u^ — 



-uf <"' = 0.3, u. 



A{0) 



B(0) 



= 



for the thermal case. The common parameters 



used in (b) are n^ = n'' = 1, m^ = 100, m^ = 1, t^^ = r^^ = r^^ = r^-^ = 1, u^^""^ = -nf (»' = 0.3, u^<°' = wf '°> = 0. 
Additionally, T = 1 for the isothermal case; T-*'"' = 10, T^*"' = 0.1 for one and T-**"' = 0.1, T^'"' = 10 for the other thermal 
cases. The second superscript "(0)" denotes the initial values. Solid lines in the figure possess corresponding theoretical slopes. 
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FIG. 2: Horizontal velocity profiles along a vertical line for the two components at time t = 2.9. The symbols denote simulation 
results from LBMs A and B. The solid line shows the analytical result. Parameters used in the simulations are m = m^ = 1, 
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= n^ = 1, 7 = 0.001. Parameters used in Eq. Ill38t are t) = 0.05, p = 1. 
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FIG. 3: LBM simulation on a diffusion process. A set of density profiles of component A are shown in (a) and the corresponding 
diffusion velocities are shown in (b). The initial density profiles of the two components follow Ea. p.'i9^ with Ap = 0.1. All the 
relaxation times are taken to be 2 x 10~^. The integration steps are Ax — Ay = 10~^ and At = 10~^. The density profiles 
from the simulation (symbols) agree well with the theoretical ones (lines). The simulation tool for this figure is LBM A. 



